****************************************************************************
**		Authors: 	Or Tuttnauer & Liran Harsgor
**		Purpose: 	Produce Figure A2 in Kedar, Harsgor & Tuttnauer (JOP)  
**		input:		KHT_countrylevel.dta
*****************************************************************************

log using "figure a2.log", replace

set more off
use "KHT_countrylevel.dta", clear

*** interaction terms
gen x1x2	=enpv * lnmeddm  	/* interaction lnmeddm and enpv*/
gen x1x3	=enpv 			* sr
gen x2x3	=		lnmeddm * sr
gen x1x2x3	=enpv * lnmeddm * sr

*** BASELINE MODEL

regress enps enpv lnmeddm x1x2 smd mmm fused upper

* Save coefficients and V-CV matrix *	
	matrix b=e(b)
	matrix V=e(V)
	scalar b1=b[1,1]  		/* first coefficient (enpv) */
	scalar b3=b[1,3]		/* interaction's coefficient (vlavgleg) */
	scalar varb1=V[1,1]		/* variance of b1 */
	scalar varb3=V[3,3]		/* variance of b3 */
	scalar covb1b3=V[1,3]	/* covariance of b1 and b3 */

	generate MVZ= ((_n-1)/14)
* Marginal effect of enpv on enps, and SE
	gen pred_meddm1=b1+b3*lnmeddm
	gen conbx= b1+b3*MVZ
	gen consx= sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ)
* Confidence interval
	gen a= 1.96*consx
	gen upperx= conbx+a
	gen lowerx= conbx-a

** Generate rugplot for scatterplot
	gen disp=0
	gen dispmeddm=lnmeddm+(runiform()-.5)/10
	gen pipe="|"

graph twoway 	 scatter disp dispmeddm if lnmeddm<3.1 & smd==0, ms(none) mlabel(pipe) mlabposition(0) mcolor(black) ///
			||	 line conbx  MVZ if MVZ<3.7, clwidth(medium) clpattern(solid) clcolor(gs10)  ///
			||   line upperx MV if MVZ<3.7, clpattern(solid) clwidth(thin) clcolor(gs10)  ///
			||   line lowerx MV if MVZ<3.7, clpattern(solid) clwidth(thin) clcolor(gs10)  /// 
            , legend(off) graphregion(fcolor(white) lcolor(white)) ///
            xtitle("Median district magnitude (logged)", size(3)) ///
            ytitle("Conversion of ENPV to ENPS", size(3))  ///
            scheme(s2mono) yscale(range(0.2 1.1)) xscale(range(0 3.7)) ylabel(0(.2)1) ///
			xlabel(0(1)3)
graph save triple_baseline.gph, replace
drop conbx consx a upperx lowerx MVZ


*** TRIPLE INTERACTION MODEL ***

*** 10% and 90% of seat ratio
sum sr, detail
scalar sr10=r(p10)
scalar sr90=r(p90)

set obs 2035000
generate sr_f = (int((_n-1)/3700))/485
gen MVZ = (mod(_n-1,3700))/1000
reg enps enpv lnmeddm sr x1x2 x1x3 x2x3 x1x2x3 smd mmm fused upper
matrix coeffs_2=e(b)

* Save coefficients and V-CV matrix *	
	matrix b=e(b)
	matrix V=e(V)
	scalar b1=b[1,1]  		/* first coefficient (enpv) */
	scalar b4=b[1,4]		/* first interaction coefficient */
	scalar b5=b[1,5]		/* second interaction coefficient */
	scalar b7=b[1,7]		/* triple interaction coefficient */
	scalar varb1=V[1,1]		/* variance of b1 */
	scalar varb4=V[4,4]		/* variance of b4 */
	scalar varb5=V[5,5]		/* variance of b5 */	
	scalar varb7=V[7,7]		/* variance of b7 */
	scalar covb1b4=V[1,4]	/* cov(b1 b4) */
	scalar covb1b5=V[1,5]	/* cov(b1 b5) */
	scalar covb4b5=V[4,5]	/* cov(b4 b5) */
	scalar covb1b7=V[1,7]	/* cov(b1 b7) */
	scalar covb4b7=V[4,7]	/* cov(b4 b7) */
	scalar covb5b7=V[5,7]	/* cov(b5 b7) */
	
** Compute predicted marginal effect of enpv for 10th and 90th percentile
** of SR, as well as the standard error of each. 
	gen conbx10= b1+b4*MVZ+b5*sr10+b7*MVZ*sr10
	gen consx10= sqrt(varb1  +varb4*(MVZ^2)  +(sr10^2)*varb5  +((sr10*MVZ)^2)*varb7 ///
	+2*covb1b4*MVZ  +2*sr10*covb1b5  +2*MVZ*sr10*covb1b7  +2*MVZ*sr10*covb4b5 ///
	+2*(MVZ^2)*sr10*covb4b7  +2*MVZ*(sr10^2)*covb5b7)
	gen conbx90= b1+b4*MVZ+b5*sr90+b7*MVZ*sr90
	gen consx90= sqrt(varb1  +varb4*(MVZ^2)  +(sr90^2)*varb5  +((sr90*MVZ)^2)*varb7 ///
	+2*covb1b4*MVZ  +2*sr90*covb1b5  +2*MVZ*sr90*covb1b7  +2*MVZ*sr90*covb4b5 ///
	+2*(MVZ^2)*sr90*covb4b7  +2*MVZ*(sr90^2)*covb5b7)

* Confidence intervals
	gen a10= 1.96*consx10
	gen upperx10= conbx10+a10
	gen lowerx10= conbx10-a10
	
	gen a90= 1.96*consx90
	gen upperx90= conbx90+a90
	gen lowerx90= conbx90-a90
	
* Generate a line at 0	
	gen yline=0
	
* Geerate graph
	graph twoway scatter disp dispmeddm if lnmeddm<3.7 in f/3700, ///
					ms(none) mlabel(pipe) mlabposition(0) mcolor(black) ///
		||	 scatter sr_f MVZ if sr_f>lowerx10 & sr_f<upperx10, mcolor(gs14) msize(vsmall) msymbol(o) ///
		||	 scatter sr_f MVZ if sr_f<upperx90 & sr_f>lowerx90, mcolor(gs11) msize(vsmall) msymbol(o) ///
		||	 scatter sr_f MVZ if sr_f<upperx10 & sr_f>lowerx90 & sr_f<upperx90, mcolor(gs9) msize(vsmall) msymbol(o) ///
		||	 line conbx10 MVZ if MVZ<3.7 in f/3700, clpattern(longdash) clwidth(medium) clcolor(black) ///
        ||   line conbx90 MVZ if MVZ<3.7 & conbx90<1.11 in f/3700, clpattern(shortdash) clwidth(medium) clcolor(black) ///
			, legend(off) graphregion(fcolor(white) lcolor(white))  ///
            xtitle("Median district magnitude (logged)", size(3)) ///
            scheme(s2mono) yscale(range(0.2 1.1)) xscale(range(0 3.7)) ylabel(0(.2)1) ///
			xlabel(0(1)3)

	graph save "triple interaction 2d.gph", replace

graph combine triple_baseline.gph "triple interaction 2d.gph", imargin(r=5 l=0) ///
		scheme(s1mono) xcommon
graph save "figure a2.gph", replace

log close
